from __future__ import division
import numpy as np
#from numpy import multiply as nm
from ShallowDonor import Shallow_Donor
from Bulk import Silicon
from Constant_Physics import Physics
from Tools import *
from MRTools import *
from MRTools import Lorentzian as lorentz
import matplotlib.pyplot as plt
import time
from scipy.integrate import odeint
import os
import src_fortran.hamiltonians as hamiltonians
import src_fortran.kronecker as kronecker
import src_fortran.ratesmatrix as ratesmatrix

from termcolor import colored
from scipy.optimize import curve_fit
from inspect import stack #To "automatize the clk functions"

#To Do list
#    - B sweep, with good linewidth
#    - fitting parameters to a B/f sweep
#    - export (.csv/.ps/...)

# Time consumption (ms)
# Parametrized_SDRmodel:
#    - EigDpRz:    6
#    - EigDpRm:    2.2
#    - EigDzRz:    8.24
#        + linalg: 0.68
#        + inv:    0.63
#        + sk22:   5.4


#New version
##########################################
##########################################
class SDR_parameters:
    #SDR model parameters
    PE = 1e8
    Cp = 1e6
    Cap = 1e6 
    Rp = 4e3
    Rap = 2e5
    
    #Modulation
    modulation_type="field" #field/AM/FM
    modulation_frequency=100000 #Hz
    modulation_phase=0 #rad
    modulation_shape="sin" #sin/square
    modulation_amplitude=0.001 #T/ percent/ MHz

    #Signal properties        
    linewidth_MHz=1
    linewidth_mT=1
    
    #Experimental conditions
    temperature=16#K 

##########################################
##########################################
class Simulation_parameters:
    time_range={}
    time_range['without resonance']=5e-6
    time_range['without resonance steady state']=3e-3
    time_range['with resonance']=300e-6
    time_range['with resonance steady state']=400e-6
    time_intervale={}
    time_intervale['without resonance']=1e-9
    time_intervale['without resonance steady state']=1e-6
    time_intervale['with resonance']=1e-6
    time_intervale['with resonance steady state']=1e-6
    do_adjust_time=False

##########################################
##########################################
class Bfield(Physics,SDR_parameters):
    
    def __init__(self,clk,B0=0.1,B0xyz=np.array([0, 0, 1]),modulated=True,fortran_compiled=True):
        self.modulated=modulated
        self.clk=clk
        self.fortran_compiled=fortran_compiled
        self.B0=B0
        self.B0xyz=B0xyz
    
    def field_vec(self,t=0):
        """
        Input:
            -t: the time
        Output:
            -field vector
        """
        res=np.array([0,0,0])
        if self.fortran_compiled==True:
            if self.modulated:
                res= hamiltonians.bfield(self.B0*Physics.FieldUnit,
                                 self.B0xyz,
                                 SDR_parameters.modulation_frequency,
                                 SDR_parameters.modulation_phase,
                                 SDR_parameters.modulation_amplitude*Physics.FieldUnit,
                                 t)
            else:
                res= hamiltonians.bfield(self.B0*Physics.FieldUnit,
                                 self.B0xyz,
                                 1,
                                 0,
                                 0,
                                 0)
        else:
            if self.modulated:
                res = (self.B0 + SDR_parameters.modulation_amplitude*Physics.FieldUnit* np.sin(2*np.pi*t*SDR_parameters.modulation_frequency+SDR_parameters.modulation_phase)) / np.linalg.norm(self.B0xyz) * self.B0xyz
            else:
                res = (self.B0 ) / np.linalg.norm(self.B0xyz) * self.B0xyz
            
        return res
        
    def field_norm(self,t=0):
        """
        Input:
            -t: the time
        Output:
            -field norm
        """
        self.clk.TIC(self.__class__.__name__+":"+stack()[0][3],"L")
        if self.modulated:
            return (self.B0 + SDR_parameters.modulation_amplitude* np.sin(2*np.pi*t*SDR_parameters.modulation_frequency+SDR_parameters.modulation_phase))
        else:
            return self.B0
        self.clk.TOC(self.__class__.__name__+":"+stack()[0][3])

class ResonanceFrequency(Physics,SDR_parameters):
    
    def __init__(self,clk,frequency,power,modulated=False):
        """
        Input:
            -f0:
            -mod_freq:
            -mod_phase:
            -mod_form:
            -mod_amp:
        """
        self.modulated=modulated
        self.clk=clk
        self.frequency=frequency
        self.power=power
        
    def rf(self,t):
        """
        Input:
            -t: the time
        Output:
            -frequency
        """
        self.clk.TIC(self.__class__.__name__+":"+stack()[0][3],"L")
        if self.modulated == False:
            return self.frequency,self.power
        elif self.modulated=="FM":
            return (self.frequency + SDR_parameters.modulation_amp* np.sin(2*np.pi*t*SDR_parameters.modulation_frequency+SDR_parameters.modulation_phase)),self.power
        elif self.modulated=="AM":
            return self.frequency,(1-SDR_parameters.modulation_amp* np.sin(2*np.pi*t*SDR_parameters.modulation_frequency+SDR_parameters.modulation_phase)/2)*self.power
        self.clk.TOC(self.__class__.__name__+":"+stack()[0][3])

##########################################
##########################################
class Lockin_Detection():
    def __init__(self,clk,fundamental_f=100000, n_harmonics=1):
        """
        Inputs:
            -fundamental_f: in Hz
            -n_harmonics: integer, 1-> 2*f0, 2-> 3*f0
        """
        self.f0=fundamental_f
        self.harmonics=np.arange(1,n_harmonics+2,1)      
        self.signals_X={}
        self.signals_Y={}  
        self.signals_R={}
        self.signals_phi={}
        self.clk=clk
        
    def detect_XY(self,time_ax,signal):
        """
        Input:
            -time_ax: the corresponding time axis, np.array([n])
            -signal: np.array([n])
        """
        for harm in self.harmonics:
            f_det=harm*self.f0
            phi0 = np.sum(signal*np.sin(2*np.pi*f_det*time_ax))
            phi90 = np.sum(signal*np.sin(2*np.pi*f_det*time_ax+np.pi/2))
            self.signals_X[harm]=phi0
            self.signals_Y[harm]=phi90

    def calculate_Rphi(self):
        for harm in self.harmonics:
            self.signals_R[harm]=np.sqrt(self.signals_X[harm]**2+self.signals_Y[harm]**2)
            # arctan2(y,x)
            self.signals_phi[harm]=np.arctan2(self.signals_Y[harm],self.signals_X[harm])

    def lockin_detect(self,signal_size,intervale,signal,harmonic=2):
        def find_x(x_axis,x):
            value=abs(x_axis[0]-x)
            index=0
            for i in range(x_axis.size):
                if abs(x_axis[i]-x)<value:
                    index=i
                    value=abs(x_axis[i]-x)
            return index
        
        freq_fft=np.fft.fftfreq(signal_size,intervale)
        conductivity_fft=np.fft.fft(signal)
        x=find_x(freq_fft,SDR_parameters.modulation_frequency*harmonic)
        return freq_fft,conductivity_fft,np.abs(conductivity_fft[x]),np.angle(conductivity_fft[x])

##########################################
##########################################
class Data():
    def __init__(self,file_name="report"):
        self.file_name=file_name
        self.conductivity={}
        self.populations={}
        self.time={}
        self.current_population=0
        self.current_conductivity=0
        self.current_time=0
        
        self.time_constant={} #Results from the exponential fitting to determine whether the steady state is reached
        self.fitted_transient={}
        self.integrated_signal_lockin_R=0
        self.integrated_signal_lockin_phi=0
        
        self.operation_order={}
        self.current_operation_id=-1

        self.detected_lockin_signal=[]
        self.detected_lockin_signal_phase=[]
    
        self.delimiter=" "
        self.newline='\n'
        self.home_path=os.getcwd()+os.sep#Current path to this program
        self.data_path=self.home_path+os.sep+"data"+os.sep#For saving later on data
        self.extension=".csv"
        
    
    def Save_lockin(self,name,phase=False):
        f = open(self.data_path+name+self.extension, 'w')
        for line in self.detected_lockin_signal:
            s=""
            s=s+str(line[0])+self.delimiter+str(line[1])+self.newline
            f.write(s)
        f.close()
        if phase:
            f = open(self.data_path+name+"_phase"+self.extension, 'w')
            for line in self.detected_lockin_signal_phase:
                s=""
                s=s+str(line[0])+self.delimiter+str(line[1])+self.newline
                f.write(s)
            f.close()

        
    def Export_conductivity(self,file_name):
        print "Saving the time domain simu as: "+ file_name+ "_td" + self.extension
        print "Saving the frequency domain simu as: "+ file_name+ "_fd" + self.extension
        print "In: "+self.data_path
        
        #time domain
        f_td = open(self.data_path+file_name+ "_td"+self.extension, 'w')
        
        #time domain
        f_fd = open(self.data_path+file_name+ "_fd"+self.extension, 'w')
        
        for i in range(0,self.current_operation_id+1):
            label=self.operation_order[i]
            if "fft" in label:
                for k in range(len(self.time[label])):
                    line=str(self.time[label][k]) + self.delimiter + str(self.conductivity[label][k].real) + self.delimiter + str(self.conductivity[label][k].imag) + self.newline
                    f_fd.write(line)
            else:
                for k in range(len(self.time[label])):
                    line=str(self.time[label][k]) + self.delimiter + str(self.conductivity[label][k]) + self.newline
                    f_td.write(line)
                
        f_td.close()
        f_fd.close()

        

        
##########################################
##########################################
class Transient(Simulation_parameters):

    def __init__(self,sdr_model,data,clk,compiled_functions=True):
        self.clk=clk
        self.data=data
        self.sdr_model=sdr_model
        self.fortran_compiled_functions=compiled_functions
        print Simulation_parameters.time_intervale["with resonance"]

    def is_steady_state_reached(self):
        """
        Description:
            Print whether the steady state is reached or not
        Todo:
            -Redefine the steady state
        Input:
            -odeint_sol: solution of the system
            -time_scale=t.size
        Output:
            -None
        Modification of properties:
            -None
        Dependence:
            -Properties:
                -self.conductivity
                -self.steadystate_margin
                -self.is_steady_state_reached_print
        """
        self.CT.TIC("is_steady_state_reached")
        if self.conductivity==[]:
            print "\nPlease compute the conductivity first."
        else:
            if np.abs(self.conductivity[-1]-self.conductivity[-2])<=self.steadystate_margin:
                if self.is_steady_state_reached_print==True:
                    print "\nSteady state reached: ds/dt= ",np.abs(self.conductivity[-1]-self.conductivity[-2])
            else:
                if self.is_steady_state_reached_print==True:
                    print "\nSteady state not reached yet: ds/dt= ",np.abs(self.conductivity[-1]-self.conductivity[-2])
                self.is_steady_state_reached_flag=False
        self.CT.TOC("is_steady_state_reached")

    def Plot_transient(self,transient_type="",modulation_plot="points"):
        """
        Input:
            - transient_type:    "conductivity": plot the conductivity transient
        """
        
        # Determine the number of plot to show
        # ie: number of transients
        n_row=2
        n_col=2        
        
        #Initialize fig
        fig = plt.figure()
        
        #Plot Light transient top left panel
        id_transient=0
        id_ss=id_transient+1 
        label_transient = self.data.operation_order[id_transient]
        label_ss = self.data.operation_order[id_ss]
        ax=fig.add_subplot(n_row,n_col,1)
        #
        inter_time=np.zeros([self.data.time[label_transient].size+self.data.time[label_ss].size+1])
        inter_time[0:self.data.time[label_transient].size]=self.data.time[label_transient]
        inter_time[self.data.time[label_transient].size:-1]=self.data.time[label_ss]
        inter_conductivity=np.zeros([self.data.conductivity[label_transient].size+self.data.conductivity[label_ss].size+1])
        inter_conductivity[0:self.data.conductivity[label_transient].size]=self.data.conductivity[label_transient]
        inter_conductivity[self.data.conductivity[label_transient].size:-1]=self.data.conductivity[label_ss]
        #
        ax.plot(inter_time[0:-2],inter_conductivity[0:-2],'b')
        ax.set_title(label_transient+' + '+label_ss)
        ax.set_ylabel('SDR(a.u.)', color='b')
        for tl in ax.get_yticklabels():
            tl.set_color('b')
            
        #Plot mw transient top right panel
        id=2
        label = self.data.operation_order[id]
        ax=fig.add_subplot(n_row,n_col,2)
        ax.plot(self.data.time[label],self.data.conductivity[label],'b')
        ax.set_title(label)
        ax.set_ylabel('SDR(a.u.)', color='b')
        for tl in ax.get_yticklabels():
            tl.set_color('b')
            
        #Plot mw transient + steady state ++ field modulation bottom left panel
        id_transient=2
        id_ss=id_transient+1 
        label_transient = self.data.operation_order[id_transient]
        label_ss = self.data.operation_order[id_ss]
        ax=fig.add_subplot(n_row,n_col,3)
        #
        inter_time=np.zeros([self.data.time[label_transient].size+self.data.time[label_ss].size+1])
        inter_time[0:self.data.time[label_transient].size]=self.data.time[label_transient]
        inter_time[self.data.time[label_transient].size:-1]=self.data.time[label_ss]
        inter_conductivity=np.zeros([self.data.conductivity[label_transient].size+self.data.conductivity[label_ss].size+1])
        inter_conductivity[0:self.data.conductivity[label_transient].size]=self.data.conductivity[label_transient]
        inter_conductivity[self.data.conductivity[label_transient].size:-1]=self.data.conductivity[label_ss]
        #
        ax.plot(inter_time[0:-2],inter_conductivity[0:-2],'b')
        ax.set_title(label_transient+' + '+label_ss)
        ax.set_ylabel('SDR(a.u.)', color='b')
        for tl in ax.get_yticklabels():
            tl.set_color('b')

        self.fieldrecord=np.zeros([inter_time.size])
        for k in range(inter_time.size):
            self.fieldrecord[k]=self.sdr_model.bfield.field_norm(inter_time[k])
        ax2 = ax.twinx()
        if modulation_plot=="line":
            ax2.plot(inter_time[0:-2], self.fieldrecord[0:-2], '-r')
        else:
            ax2.plot(inter_time[0:-2], self.fieldrecord[0:-2], '.r')
        ax2.set_ylabel('Field(T)', color='r')
        for tl in ax2.get_yticklabels():
            tl.set_color('r')
            
        #Plot fft of the steady state
        id=4
        label = self.data.operation_order[id]
        ax3=fig.add_subplot(n_row,n_col,4)
        ax3.plot(self.data.time[label],self.data.conductivity[label],'b')
        ax3.set_ylabel('fft', color='b')
        ax3.set_xlim(-int(SDR_parameters.modulation_frequency/2),int(SDR_parameters.modulation_frequency*2.5))
        ax3.set_ylim(bottom=0,top=400)
        
        plt.show()
                
    def Temporal_evolution(self,label,excitation_list=[]):
        self.data.time[label]=Arange(self.data.current_time,
                                     self.data.current_time+Simulation_parameters.time_range[label],
                                     Simulation_parameters.time_intervale[label])
        self.Temporal_evolution_processing(label,
                                           do_fit=Simulation_parameters.do_adjust_time,
                                           excitation_list=excitation_list)
        if Simulation_parameters.do_adjust_time:
            if Simulation_parameters.time_range[label]<self.data.time_constant[label]*Simulation_parameters.transient_std_length+Simulation_parameters.integration_time:
                # Not enough points to integrate
                label2=label+" remaining"
                #Duplicate parameters
                Simulation_parameters.time_intervale[label2]=Simulation_parameters.time_intervale[label]
                Simulation_parameters.time_range[label2]=Simulation_parameters.time_range[label]
                
                remaining_time=self.data.time_constant[label]*Simulation_parameters.transient_std_length+Simulation_parameters.integration_time-Simulation_parameters.time_range[label]
                self.data.time[label2]=Arange(self.data.current_time,
                                             self.data.current_time+remaining_time,
                                             Simulation_parameters.time_intervale[label])
                self.Temporal_evolution_processing(label2,do_fit=False,excitation_list=excitation_list)
                self.data.time[label]=np.concatenate([self.data.time[label],self.data.time[label2]],axis=0)
                self.data.populations[label]=np.concatenate([self.data.populations[label],self.data.populations[label2]],axis=0)
                self.data.conductivity[label]=np.concatenate([self.data.conductivity[label],self.data.conductivity[label2]],axis=0)
        
            def find_t_index(t,label):
                #t=a * t_index + b
                a=Simulation_parameters.time_intervale[label]
                b=self.data.time[label][0]
                return int((t-b)/a)
            label_transient=label+" transient"
            #Duplicate parameters
            Simulation_parameters.time_intervale[label_transient]=Simulation_parameters.time_intervale[label]
            Simulation_parameters.time_range[label_transient]=Simulation_parameters.time_range[label]
            label_ss=label+" steady state"
            #Duplicate parameters
            Simulation_parameters.time_intervale[label_ss]=Simulation_parameters.time_intervale[label]
            Simulation_parameters.time_range[label_ss]=Simulation_parameters.time_range[label]
            #
            t1_index=0
            t1=self.data.time[label][t1_index]
            
            t2=self.data.time_constant[label]*Simulation_parameters.transient_std_length+t1
            t2_index=find_t_index(t2,label)
            
            t3=Simulation_parameters.integration_time+t2
            t3_index=find_t_index(t3,label)
            
    
            self.data.current_operation_id+=1          
            self.data.time[label_transient]=self.data.time[label][t1_index:t2_index]
            self.data.populations[label_transient]=self.data.populations[label][t1_index:t2_index]
            self.data.conductivity[label_transient]=self.data.conductivity[label][t1_index:t2_index]
            self.data.operation_order[self.data.current_operation_id]=label_transient
            
            self.data.current_operation_id+=1  
            self.data.time[label_ss]=self.data.time[label][t2_index+1:t3_index]
            self.data.populations[label_ss]=self.data.populations[label][t2_index+1:t3_index]
            self.data.conductivity[label_ss]=self.data.conductivity[label][t2_index+1:t3_index]
            self.data.operation_order[self.data.current_operation_id]=label_ss
            
        else:
            label_transient=label+" transient"
            #Duplicate parameters
            Simulation_parameters.time_intervale[label_transient]=Simulation_parameters.time_intervale[label]
            Simulation_parameters.time_range[label_transient]=Simulation_parameters.time_range[label]
            
            self.data.current_operation_id+=1          
            self.data.time[label_transient]=self.data.time[label]
            self.data.populations[label_transient]=self.data.populations[label]
            self.data.conductivity[label_transient]=self.data.conductivity[label]
            self.data.operation_order[self.data.current_operation_id]=label_transient
        
    def Temporal_evolution_processing(self,label,do_fit=False,excitation_list=[]):
        self.clk.TIC(self.__class__.__name__+":"+stack()[0][3]+":"+label,"H")
        #Print information
        print colored("\nTemporal evolution "+label+":",'blue')
    
               
        # defining the function dy/dt = f(y,t)

            
        def func(y,t0):
            if self.data.current_time!=t0:
                
                self.clk.TIC("Do Eig","L")
                self.sdr_model.Do_Eig(t0)
                self.clk.TOC("Do Eig")
                
                self.clk.TIC("fr.nonresmat","L")
                #self.RatesMatrix=self.sdr_model.NonResonantMatrix()
                self.sdr_model.NonResonantMatrix()
                self.clk.TOC("fr.nonresmat")
                
                self.clk.TIC("fr.resmat","L")
                #self.RatesMatrix+=self.sdr_model.ResonantMatrix(excitation_list,self.data.current_time)
                self.sdr_model.ResonantMatrix(excitation_list,self.data.current_time)
                self.clk.TOC("fr.resmat")
                    
            self.data.current_time=t0
            Print_progress(t0,self.data.time[label][-1],start=self.data.time[label][0])
            #print np.sum(y)
            return np.dot(self.sdr_model.RatesMatrix,y)

        def Dfunc(y,t0):
            
            if self.data.current_time!=t0:
                self.clk.TIC("Do Eig","L")
                self.sdr_model.Do_Eig(t0)
                self.clk.TOC("Do Eig")
                
                self.clk.TIC("fr.nonresmat","L")
                #self.RatesMatrix=self.sdr_model.NonResonantMatrix()
                self.sdr_model.NonResonantMatrix()
                self.clk.TOC("fr.nonresmat")
                
                self.clk.TIC("fr.resmat","L")
                #self.RatesMatrix+=self.sdr_model.ResonantMatrix(excitation_list,self.data.current_time)
                self.sdr_model.ResonantMatrix(excitation_list,self.data.current_time)
                self.clk.TOC("fr.resmat")
            
            self.data.current_time=t0
        
            return self.sdr_model.RatesMatrix
        
        # performing the ODEINT
        y=odeint(func,
                 self.data.current_population,
                 self.data.time[label],
                 #hmax=Simulation_parameters.time_intervale[label],
                 #hmax=SDR_parameters.modulation_frequency/10,
                 #hmin=SDR_parameters.modulation_frequency/10000,
                 Dfun=Dfunc,
                 mxstep=2000)
  
        #taking the last solution step (last line) and save it as steadystate, and storing the corresponding conductivity
        self.data.populations[label]=y
        self.data.conductivity[label]=self.sdr_model.Compute_conductivity(y,self.data.time[label].size)
        self.data.current_population = y[-1,:]
        self.data.current_conductivity=self.data.conductivity[label][-1]
        
        #Fitting with an exponential and save to data
        if do_fit==True:
            self.fit_transient(label)
        
        self.clk.TOC(self.__class__.__name__+":"+stack()[0][3]+":"+label)
        
    def fit_transient(self,label):
        """
        We suppose ft to be exp
        We want to know "tau"
        """
        self.clk.TIC(self.__class__.__name__+":"+stack()[0][3],"M")
        t=self.data.time[label]-self.data.time[label][0]
        ft=self.data.conductivity[label]-self.data.conductivity[label][0]
        def func(t,A,tau):
            return A*(1-np.exp(-t/tau))
        popt,pcov=curve_fit(func,t,ft)
        self.data.time_constant[label]=popt[1]
        self.data.fitted_transient[label]=popt[0]*(1-np.exp(-t/popt[1]))
        self.clk.TOC(self.__class__.__name__+":"+stack()[0][3])
        return popt
    
##########################################
##########################################

class Parametrized_SDR_model(SDR_parameters):
    
    def __init__(self, bulk, donor,bfield,clk,excitation_dic={},fortran_compiled=True):
        """
        specie: define tje shallow donor specie ("P", "As", "Sb5/2", "Sb7/2", "Bi")
        D: Dissociation rate
        G: Generation rate
        C: Capture rate
        Rp: Recombination rate of parallel pairs
        Rap: Recombination rate of antiparallel pairs
        W: carrier generation due to illumination
        temperature: in Kelvin
        average_dist: average separation between one donor and one PRC, in 'nm'
        """
        self.clk=clk
        self.fortran_compiled_functions=fortran_compiled
        
        self.bfield=bfield
        self.excitation_dic=excitation_dic
        self.donor = donor
        self.bulk = bulk

        #-----------------
        self.dim_een = int(self.bulk.dim_PRC * self.donor.dim_e * self.donor.dim_n) #number of states with the donor ionized
        self.full_dim = self.dim_een + self.dim_een + self.donor.dim_n

        self.Mbj = np.zeros([self.dim_een, self.dim_een])
        self.Mjb = np.zeros([self.dim_een, self.dim_een])
        self.Mbj_el_wise_square = np.zeros([self.dim_een, self.dim_een])
        #-----------------
        
        # Rate matrix: column names
            #[0...c1e, c2i=c1e...c2e, c3i=c2e...c3e]
            #[                                                ]
            #[                                                ]
            #[                                                ]
            #[                                                ]
        self.c1i=0
        self.c1e=self.dim_een
        self.c2i=self.c1e
        self.c2e=self.dim_een + self.dim_een
        self.c3i=self.c2e
        self.c3e=self.full_dim
        
        
        self.dipolar_transitions=[]    # [[B0,li-lj,EPR,SDR],...]

        self.population_away_from_resonance = np.zeros(self.dim_een)
        # Resonant b term
        self.population_at_resonance = np.zeros(self.dim_een)
        self.RatesMatrix=np.zeros([self.dim_een+self.dim_een+self.donor.dim_n,self.dim_een+self.dim_een+self.donor.dim_n])


        #Construction of the eigenvectors lists and hamiltonians matrix
        self.eigvalDpRz = np.zeros([self.dim_een])
        self.eigvalDzRz = np.zeros([self.dim_een])
        self.hamiltonianDzRz = np.zeros([self.dim_een, self.dim_een])
        self.eigvalDpRm = np.zeros([self.donor.dim_n])
                
    def EigDpRz(self,t):
        """
        Description:
            Define the Hamiltonian in the DpRz state
        Todo:
            -None
        Input:
            -B: magnetic field
        Output:
            -self.eigvalDpRz: matrix of the Hamiltonian in the case D+ R0
        Modification of properties:
            -self.hamiltonianDpRz
        Dependence:
            -Properties:
                -self.dim_een
                -self.B
            -Class:
                -self.donor
                -self.bulk
        Comments:
            -relaxation supposed to be faster than others processes
        """
        #nuclear Zeeman and PRC
        #we write only the eigenvalues as it is already diagonalized Assuming B=Bz (here, avoid a diagonalization)
        self.clk.TIC(self.__class__.__name__+":"+stack()[0][3],"L")
        #self.hamiltonianDpRz=np.kron(np.kron(self.bulk.PRC_Hamiltonian_Ze(self.bfield.field_vec(t)),self.donor.Sd),self.donor.Id)
        #self.hamiltonianDpRz+=np.kron(np.kron(self.bulk.Sd,self.donor.D_Hamiltonian_Ze(self.bfield.field_vec(t))),self.donor.Id)
        #self.hamiltonianDpRz+=np.kron(np.kron(self.bulk.Sd,self.donor.Sd),self.donor.D_Hamiltonian_Zn(self.bfield.field_vec(t)))
        bfieldvec=self.bfield.field_vec(t)
        o=np.zeros([self.full_dim,self.full_dim])
        self.hamiltonianDpRz=kronecker.sskk22(o,
                                       self.bulk.dim_PRC,
                                       self.donor.dim_e,
                                       self.donor.dim_n,
                                       self.bulk.PRC_Hamiltonian_Ze(bfieldvec),
                                       self.donor.D_Hamiltonian_Ze(bfieldvec),
                                       self.donor.D_Hamiltonian_Zn(bfieldvec),
                                       self.bulk.Sd,
                                       self.donor.Sd,
                                       self.donor.Id)
        self.clk.TOC(self.__class__.__name__+":"+stack()[0][3])

    def EigDzRz(self,t):
        """
        Description:
            Define the Hamiltonian in the DzRz state
        Todo:
            -None
        Input:
            -None
        Output:
            -self.eigvalDpRz: matrix of the Hamiltonian in the case D+ R0
        Modification of properties:
            -self.hamiltonianDpRz
            -self.Mbj
            -self.Mjb
            -self.Mbj_el_wise_square
            -self.Mjb_el_wise_square
        Dependence:
            -Properties:
                -self.dim_een
                -self.B
                -self.coupling_PRC_donor: Allows B(t) or not
            -Class:
                -self.donor
                -self.bulk
        Comments:
            -relaxation supposed to be faster than others processes
        """
        #self.hamiltonianDzRz = np.kron(self.bulk.Sd,o)    #Donor term
        #self.hamiltonianDzRz += np.kron(self.bulk.PRC_Hamiltonian_Ze(self.bfield.field_vec(t)),np.kron(self.donor.Sd, self.donor.Id)) #PRC term
        b=self.bfield.field_vec(t)
        if self.fortran_compiled_functions==True:
            self.hamiltonianDzRz=kronecker.sk22( self.bulk.dim_PRC,
                                          self.donor.dim_e*self.donor.dim_n,
                                          1,1,1,0,1,1,
                                          self.bulk.PRC_Hamiltonian_Ze(b),
                                          self.donor.D_Hamiltonian_Dz(b),
                                          self.bulk.Sd,
                                          kronecker.k22(self.donor.Sd,self.donor.dim_e,self.donor.dim_e,1, self.donor.Id,self.donor.dim_n,self.donor.dim_n,1,1),
                                          1)
            
        else:
            self.hamiltonianDzRz=np.kron(self.bulk.PRC_Hamiltonian_Ze(b),np.kron(self.donor.Sd,self.donor.Id))+np.kron(self.bulk.Sd,self.donor.D_Hamiltonian_Dz(b))
        
        #because in the end, this is a real matrix
        self.hamiltonianDzRz=np.real(self.hamiltonianDzRz)
        self.eigvalDzRz, self.Mbj = np.linalg.eig(self.hamiltonianDzRz)  #Diagonalization
        self.Mjb = np.linalg.inv(self.Mbj)
        
        self.Mbj_el_wise_square=self.Mbj**2
        
    def EigDpRm(self,t):
        """
        Description:
            Define the Hamiltonian in the DzRz state
        Todo:
            -None
        Input:
            -B: magnetic field
        Output:
            -None
        Modification of properties:
            -self.hamiltonianDpRm
            -self.Mkb
            -self.Mbk
        Dependence:
            -Properties:
                -self.dim_een
                -self.B
            -Class:
                -self.donor
                -self.bulk
        Comments:
            -relaxation supposed to be faster than others processes
        """
        #nuclear Zeeman and PRC
        #we write only the eigenvalues as it is already diagonalized
        self.hamiltonianDpRm=self.donor.D_Hamiltonian_Zn(self.bfield.field_vec(t))
        
    def Do_Eig(self,t=0):
        """
        Description:
            Perform the eigensystem, fill the transitions list...
            20120911,9am: execution time around 2 ms
        Todo:
            -None
        Input:
            -B: from MRTools.Bfield(B,xyz)
        Output:
            -None
        Modification of properties:
            -self.Dipolar_Transition_Matrix
            -self.n_of_transitions
            -self.dipolar_transitions
            -self.nonzero_i_index
            -self.nonzero_j_index
        Dependence:
            -Properties:
                -self.Mjb
                -self.Mbj
                -self.Dipolar_Transition_Matrix
            -Methods:
                -self._EigDpRz()
                -self._EigDzRz()
                -self._EigDpRm()
            -class:
                -self.donor
                -self.bulk
        """
        
        #self.EigDpRz(t)
        self.EigDzRz(t)
        #self.EigDpRm(t)
        
        DP1=np.kron(np.kron(self.bulk.Sd,self.donor.Sx),self.donor.Id)    #Dipolar transition
        self.Dipolar_Transition_Matrix=np.dot(self.Mjb,np.dot(DP1,self.Mbj))   #Change the basis
        self.n_of_transitions=np.count_nonzero(self.Dipolar_Transition_Matrix)
        self.dipolar_transitions=np.zeros([self.n_of_transitions,4])    # [[B0,li-lj,EPR],...]
        self.nonzero_i_index,self.nonzero_j_index=np.nonzero(self.Dipolar_Transition_Matrix)

    def NonResonantMatrix(self):
                
        if self.fortran_compiled_functions==True:
            self.RatesMatrix= ratesmatrix.nonresmat(self.Mbj_el_wise_square,
                                          self.dim_een+self.dim_een+self.donor.dim_n,
                                          self.dim_een,
                                          self.bulk.dim_PRC,
                                          self.donor.dim_e,
                                          self.donor.dim_n,
                                          SDR_parameters.Cp,
                                          SDR_parameters.Cap,
                                          SDR_parameters.Rp,
                                          SDR_parameters.Rap,
                                          SDR_parameters.PE)
        else:
                # Ci: assuming there is no coupling at this step (even if we write the transformation matrix)
            Calpha = (SDR_parameters.Cp - SDR_parameters.Cap) / 2
            Cbeta = (SDR_parameters.Cp + SDR_parameters.Cap) / 2
            Ci_python = np.zeros([self.dim_een, self.dim_een])
            Ci_python = Calpha * np.kron(np.kron(self.bulk.Sz, self.donor.Sz), self.donor.Id) + Cbeta * np.kron(np.kron(self.bulk.Sd, self.donor.Sd), self.donor.Id)
            
                # Cj
            Cj_python = np.zeros([self.dim_een, self.dim_een])
            v=np.kron(np.array([SDR_parameters.Cp,SDR_parameters.Cap,SDR_parameters.Cap,SDR_parameters.Cp]),np.ones([self.donor.dim_n]))
            for k in range(self.dim_een):
                Cj_python[k,:]=(v*self.Mbj_el_wise_square[:,k]).T
                #print v.size
                #print (v*self.Mbj_el_wise_square[:,k]).size
                #print Cj[k,:].size
            
                # Rj
            Rj_python = np.zeros([self.dim_een, self.dim_een])
            R2=np.array([SDR_parameters.Rp,SDR_parameters.Rap,SDR_parameters.Rap,SDR_parameters.Rp])
            R3=np.ones([self.donor.dim_n])
            R1=np.kron(R2,R3)
            for j in range(self.dim_een):
                Rj_python[j,j]=np.dot(R1.T,self.Mbj_el_wise_square[:,j])
            
                # Rk
            Rk_python = np.zeros([self.donor.dim_n, self.dim_een])
            for j in range(self.donor.dim_n):
                R4=np.zeros([self.donor.dim_n])
                R4[j]=1
                R4=np.kron(R2,R4)
                for i in range(self.dim_een):
                    Rk_python[j,i]=np.dot(R4,self.Mbj_el_wise_square[:,i])
            
            
                # W: we assume there is no polarization of the excitation light ...
            PEi_python = SDR_parameters.PE/self.dim_een*np.kron(np.ones([self.bulk.dim_PRC*self.donor.dim_e,1]),self.donor.Id)
            PEk_python = SDR_parameters.PE / self.donor.dim_n * self.donor.Id
            
            res=np.zeros([self.c3e,self.c3e])
            self.RatesMatrix[self.c1i:self.c1e, self.c1i:self.c1e] = -Ci_python
            self.RatesMatrix[self.c2i:self.c2e, self.c1i:self.c1e] = Cj_python
            self.RatesMatrix[self.c2i:self.c2e, self.c2i:self.c2e] = -Rj_python
            self.RatesMatrix[self.c3i:self.c3e, self.c2i:self.c2e] = Rk_python
            self.RatesMatrix[self.c1i:self.c1e, self.c3i:self.c3e] = PEi_python
            self.RatesMatrix[self.c3i:self.c3e, self.c3i:self.c3e] = -PEk_python
    
    def ResonantMatrix(self,excitation_list,current_time):
        res=np.zeros([self.c3e,self.c3e])
        if self.fortran_compiled_functions==True:
            for exc in excitation_list:
                f,power=self.excitation_dic[exc].rf(current_time)
                self.RatesMatrix=ratesmatrix.resmat(self.RatesMatrix,
                                                    self.Dipolar_Transition_Matrix,
                                                    f,
                                                    self.linewidth_MHz,
                                                    power,
                                                    self.eigvalDzRz,
                                                    self.dim_een+self.dim_een+self.donor.dim_n,
                                                    self.dim_een)
        else:
            
            lorentzian = lorentz(derivative=0)
                
            for i in range(self.Dipolar_Transition_Matrix[:, 0].size):
                for j in range(self.Dipolar_Transition_Matrix[0, :].size):
                    for exc in excitation_list:
                        f,power=self.excitation_dic[exc].rf(current_time)
                        center = np.abs(self.eigvalDzRz[i] - self.eigvalDzRz[j])
                        intensity = self.Dipolar_Transition_Matrix[i, j] ** 2
                        #dipole_transition_strength = lorentzian.Eval(coeffs, evalpts)
                        dipole_transition_strength = ratesmatrix.lorentzian(intensity,f,self.linewidth_MHz,center)
                                            
                        im=self.dim_een+i
                        jm=self.dim_een+j
                    
                        self.RatesMatrix[im,im]+=-power*dipole_transition_strength
                        self.RatesMatrix[jm,jm]+=-power*dipole_transition_strength
                        self.RatesMatrix[im,jm]+=power*dipole_transition_strength
                        self.RatesMatrix[jm,im]+=power*dipole_transition_strength

    def Compute_conductivity(self,odeint_sol,time_scale):
        """
        Description:
            Compute the conductivity knowing the populations
        Todo:
            -None
        Input:
            -odeint_sol: solution of the system
            -time_scale=t.size
        Output:
            -self.conductivity
        Modification of properties:
            -self.conductivity
        Dependence:
            -Properties:
                -self.c1i,self.c1e
                -self.c3i,self.c3e
            -class:
                -self.donor
                -self.bulk
        
        """
        self.conductivity=np.zeros([time_scale])
        for k in range(time_scale):
            self.conductivity[k]=(self.bulk.mue-self.bulk.muh)*np.sum(odeint_sol[k,self.c1i:self.c1e])
        return self.conductivity
          
    def Check_list(self):
        """
        Description:
            Print whether the matrixes are population conservative or not, and if the steady state is reached or not
        Todo:
            -None
        Input:
            -None
        Output:
            -None
        Modification of properties:
            -None
        Dependence:
            -Properties:
                -self.non_resonant_matrix_pop_conservative
                -self.resonant_matrix_pop_conservative
                -self.is_steady_state_reached_flag
        """
        self.clk.TIC("Check_list")
        if not self.non_resonant_matrix_pop_conservative:
            print "\nThe non resonant matrix does not conserve the population."
        
        if not self.resonant_matrix_pop_conservative:
            print "\nAt least one of the resonant matrix does not conserve the population."
        
        if not self.is_steady_state_reached_flag:
            print "\nOne of the steady states at least was not reached."
        self.clk.TOC("Check_list")
  
    def Plot_populations(self,ploted_variable="steady state population",plot_size="full"):
        """
        Input:
            - ploted_variable=      "steady state populations"
                                    "resoning state populations"
                                    "difference populations"
            - plot_size =   "full"
                            "donor"
        """
        self.clk.TIC("Plot_populations")
        if plot_size=="full":
            ind1 = range(self.full_dim)  # the x locations for the groups
        elif plot_size=="donor":
            ind1 = range(self.c2i,self.c2e)  # the x locations for the groups
            
        w = 0.4       # the width of the bars
        y1=[]
        for k in ind1:
            if ploted_variable =="steady state populations":
                y1.append(self.steadystate[k])
                title="Steady state populations"
            elif ploted_variable =="resoning state populations":
                y1.append(self.resoningstate[k])
                title="Populations at resonance"
            elif ploted_variable =="difference populations":
                y1.append(self.steadystate[k]-self.resoningstate[k])
                title="Populations at SS - at resonance"
        
        fig=plt.figure()
        ax = fig.add_subplot(1,1,1)
        
        ax.bar(ind1, y1, width=w,color='r')
        
        # add some
        plt.ylabel('Populations')
        plt.title(title)
        
        plt.show()
        self.clk.TOC("Plot_populations")

class Field_sweep():
    def __init__(self):
        #General instanciation
        self.clk=Computation_time(True)
        self.data=Data("Report")
        self.sdr_param=SDR_parameters()
        self.Bi = Shallow_Donor(self.clk,"Bi")
        self.P = Shallow_Donor(self.clk,"P")
        self.Si = Silicon(self.clk)
        self.lockin_base=Lockin_Detection(self.clk)
        self.lockin_resonance=Lockin_Detection(self.clk)
        
        
        self.rf1=ResonanceFrequency(self.clk,frequency=9076,power=1e10,modulated=False)
        self.excitation_dic={'rf1':self.rf1}
    
        
        #Specific instanciation
        self.field_list=[0.0660558,
                         0.0796147,
                         0.0986649,
                         0.12595,
                         0.164827,
                         0.217885,
                         0.285122,
                         0.363939,
                         0.450981,
                         0.543407]
        self.intensity_harm1=[]
        self.intensity_harm2=[]
        self.phase_harm1=[]
        self.phase_harm2=[]
        
    def Do_sweep(self):
        self.clk.TIC(self.__class__.__name__+":"+stack()[0][3],"H")
        for b in self.field_list:
            print '\nMagnetic field B0= ', b
            #Specific instanciation
            self.field=Bfield(self.clk,B0=b,modulated=True)        
            self.sdr_model=Parametrized_SDR_model(bulk=self.Si, donor=self.Bi,bfield=self.field, clk=self.clk)
            self.transient=Transient(sdr_model=self.sdr_model,data=self.data,clk=self.clk,excitation=self.excitation)
            
            #########################
            #Initialization
            self.data.current_population=np.zeros([self.sdr_model.full_dim])
            self.data.current_population[self.sdr_model.c3i:self.sdr_model.c3e]=1/self.sdr_model.donor.dim_n*np.ones([self.sdr_model.donor.dim_n])
            
            #########################
            #Time evolution
            self.transient.Temporal_evolution_time_adjusted('without resonance')
            self.lockin_base.detect_XY(self.data.time['without resonance transient'], self.data.conductivity['without resonance transient'])
            self.lockin_base.calculate_Rphi()
            
            self.transient.Temporal_evolution_time_adjusted('with resonance',excitation_list=['rf1'])
            self.lockin_resonance.detect_XY(self.data.time['with resonance transient'], self.data.conductivity['with resonance transient'])
            self.lockin_resonance.calculate_Rphi()
            #print lockin_resonance.signals_R[1],lockin_resonance.signals_R[2]
            
            self.intensity_harm1.append(self.lockin_resonance.signals_R[1]-self.lockin_base.signals_R[1])
            self.intensity_harm2.append(self.lockin_resonance.signals_R[2]-self.lockin_base.signals_R[2])
            self.phase_harm1.append(self.lockin_resonance.signals_phi[1])
            self.phase_harm2.append(self.lockin_resonance.signals_phi[2])
            print "\nWithout resonance: R1= ",self.intensity_harm1[-1]," R2= ",self.intensity_harm2[-1]
            print "Without resonance: phi1= ",self.phase_harm1[-1]," phi2= ",self.phase_harm2[-1]
        
        self.clk.TOC(self.__class__.__name__+":"+stack()[0][3])
        self.clk.Display(level="H") #L/M/H
    
    def Plot_sweep(self):
        n_plots=2
        
        #Initialize fig
        fig = plt.figure()
        
        #Plot harmonic 1
        ax=fig.add_subplot(2,1,1)
        ax.plot(self.field_list,self.intensity_harm1,'ob-')
        ax.set_title("1st harmonic")
        ax.set_ylabel('SDR(a.u.)', color='b')
        for tl in get_yticklabels():
            tl.set_color('b')

        ax2 = ax.twinx()
        ax2.plot(self.field_list, self.phase_harm1, 'or-')
        ax2.set_ylabel('Phi', color='r')
        for tl in ax2.get_yticklabels():
            tl.set_color('r')
        
        #Plot harmonic 1
        ax=fig.add_subplot(2,1,2)
        ax.plot(self.field_list,self.intensity_harm2,'ob-')
        ax.set_title("2nd harmonic")
        ax.set_ylabel('SDR(a.u.)', color='b')
        for tl in ax.get_yticklabels():
            tl.set_color('b')

        ax2 = ax.twinx()
        ax2.plot(self.field_list, self.phase_harm2, 'or-')
        ax2.set_ylabel('Phi', color='r')
        for tl in ax2.get_yticklabels():
            tl.set_color('r')
        
        plt.show()

def Main(transition=0,file_name="SDR_simulation"):
    #########################
    #Instanciation
    clk=Computation_time(True)
    clk.TIC("Total","H")
    data=Data("Report")
    sdr_param=SDR_parameters()
    Bi = Shallow_Donor(clk,"Bi",fortran_compiled=True)
    Si = Silicon(clk,fortran_compiled=True)
    lockin=Lockin_Detection(clk)
    
    field_list=[0.0660558,
                 0.0796147,
                 0.0986649,
                 0.12595,
                 0.164827,
                 0.217885,
                 0.285122,
                 0.363939,
                 0.450981,
                 0.543407]
    
    b=Bfield(clk,field_list[transition],modulated=True)
    rf1=ResonanceFrequency(clk,frequency=9076,power=1e10,modulated=False)
    rf2=ResonanceFrequency(clk,frequency=100,power=1e10,modulated=False)
    excitation={'rf1':rf1}
    #excitation={'rf1':rf1,'rf2':rf2}
    
    sdr_model=Parametrized_SDR_model(bulk=Si, donor=Bi,bfield=b, clk=clk,excitation_dic=excitation)
    transient=Transient(sdr_model=sdr_model,data=data,clk=clk)
    
    #########################
    #Initialization
    data.current_population=np.zeros([sdr_model.full_dim])
    data.current_population[sdr_model.c3i:sdr_model.c3e]=1/sdr_model.donor.dim_n*np.ones([sdr_model.donor.dim_n])
    
    #########################
    #Time evolution
    label='without resonance'
    transient.Temporal_evolution(label)
    label='without resonance steady state'
    transient.Temporal_evolution(label)
    label='with resonance'
    transient.Temporal_evolution(label,excitation_list=['rf1'])
    label='with resonance steady state'
    transient.Temporal_evolution(label,excitation_list=['rf1'])
    
    
    res_freq,res_signal,res_lockin,phase= lockin.lockin_detect(transient.data.conductivity[label][0:-2].size,
                                                         Simulation_parameters.time_intervale[label],
                                                         transient.data.conductivity[label][0:-2],
                                                         2)
    label='with resonance steady state fft'
    transient.data.current_operation_id+=1
    transient.data.operation_order[transient.data.current_operation_id]=label
    transient.data.time[label]=res_freq
    transient.data.conductivity[label]=res_signal
    transient.data.detected_lockin_signal.append([field_list[transition],res_lockin])
    #########################
    # Plot results and save
    transient.data.Export_conductivity(file_name)
    print transient.data.detected_lockin_signal
    transient.Plot_transient(modulation_plot="line")# conductivity,light and resonance transient, fit/field
 
    clk.TOC("Total")
    clk.Display(level="L") #L/M/H


def Main_modfreq(transition=0,file_name="SDR_simulation"):
    #########################
    #Instanciation
    clk=Computation_time(True)
    clk.TIC("Total","H")
    data=Data("Report")
    
    Bi = Shallow_Donor(clk,"Bi",fortran_compiled=True)
    Si = Silicon(clk,fortran_compiled=True)
    lockin=Lockin_Detection(clk)
    
    field_list=[0.0660558,
                 0.0796147,
                 0.0986649,
                 0.12595,
                 0.164827,
                 0.217885,
                 0.285122,
                 0.363939,
                 0.450981,
                 0.543407]
    

    
    
    for k in np.logspace(2,5,2):
        class Simulation_parameters:
            time_range={}
            time_range['without resonance']=5e-6
            time_range['without resonance steady state']=3e-3
            time_range['with resonance']=100/k #1000 us for 100kHz mod
            time_intervale={}
            time_intervale['without resonance']=1e-9
            time_intervale['without resonance steady state']=1e-6
            time_intervale['with resonance']=1/(10*k)
            do_adjust_time=False
            
        class SDR_parameters:
            #SDR model parameters
            PE = 1e8
            Cp = 1e6
            Cap = 1e6 
            Rp = 4e3
            Rap = 2e5
            
            #Modulation
            modulation_type="field" #field/AM/FM
            modulation_frequency=k #Hz
            modulation_phase=0 #rad
            modulation_shape="sin" #sin/square
            modulation_amplitude=0.001 #T/ percent/ MHz
        
            #Signal properties        
            linewidth_MHz=1
            linewidth_mT=1
            
            #Experimental conditions
            temperature=16#K 
            
        b=Bfield(clk,field_list[transition],modulated=True)
        rf1=ResonanceFrequency(clk,frequency=9076,power=1e10,modulated=False)
        excitation={'rf1':rf1}
        
        sdr_model=Parametrized_SDR_model(bulk=Si, donor=Bi,bfield=b, clk=clk,excitation_dic=excitation)
        transient=Transient(sdr_model=sdr_model,data=data,clk=clk)
            
        #########################
        #Initialization
        data.current_population=np.zeros([sdr_model.full_dim])
        data.current_population[sdr_model.c3i:sdr_model.c3e]=1/sdr_model.donor.dim_n*np.ones([sdr_model.donor.dim_n])
        
        print "Modulation frequency: ",k
        print "Time range: ",Simulation_parameters.time_intervale['with resonance']
        #########################
        #Time evolution
        label='without resonance'
        transient.Temporal_evolution(label)
        label='without resonance steady state'
        transient.Temporal_evolution(label)
        label='with resonance'
        transient.Temporal_evolution(label,excitation_list=['rf1'])
        
        
        res_freq,res_signal,res_lockin,phase= lockin.lockin_detect(transient.data.conductivity[label][0:-2].size,
                                                             Simulation_parameters.time_intervale[label],
                                                             transient.data.conductivity[label][0:-2],
                                                             2)
        label='with resonance steady state fft'
        transient.data.current_operation_id+=1
        transient.data.operation_order[transient.data.current_operation_id]=label
        transient.data.time[label]=res_freq
        transient.data.conductivity[label]=res_signal
        transient.data.detected_lockin_signal.append([field_list[transition],res_lockin])
        #########################
        # Plot results and save
        transient.data.Export_conductivity(file_name+"_"+str(k))
        
    print transient.data.detected_lockin_signal
    
    clk.TOC("Total")
    clk.Display(level="L") #L/M/H 
 
def Main_Spectrum():
    #########################
    #Instanciation
    clk=Computation_time(True)
    clk.TIC("Total","H")
    data=Data("Report")
    sdr_param=SDR_parameters()
    Bi = Shallow_Donor(clk,"Bi",fortran_compiled=True)
    Si = Silicon(clk,fortran_compiled=True)
    lockin=Lockin_Detection(clk)
    
    field_list=[0.0660558,
                 0.0796147,
                 0.0986649,
                 0.12595,
                 0.164827,
                 0.217885,
                 0.285122,
                 0.363939,
                 0.450981,
                 0.543407]
    for transition in range(10):
        print "\n\nField= ",field_list[transition]
        b=Bfield(clk,field_list[transition],modulated=True)
        rf1=ResonanceFrequency(clk,frequency=9076,power=1e10,modulated=False)
        rf2=ResonanceFrequency(clk,frequency=149.885,power=1e10,modulated=False)
        #excitation={'rf1':rf1}
        excitation={'rf1':rf1,'rf2':rf2}
        
        sdr_model=Parametrized_SDR_model(bulk=Si, donor=Bi,bfield=b, clk=clk,excitation_dic=excitation)
        transient=Transient(sdr_model=sdr_model,data=data,clk=clk)
        
        #########################
        #Initialization
        data.current_population=np.zeros([sdr_model.full_dim])
        data.current_population[sdr_model.c3i:sdr_model.c3e]=1/sdr_model.donor.dim_n*np.ones([sdr_model.donor.dim_n])
        
        #########################
        #Time evolution
        label='without resonance'
        transient.Temporal_evolution(label)
        label='with resonance'
        transient.Temporal_evolution(label,excitation_list=['rf1','rf2'])
        label='with resonance steady state'
        transient.Temporal_evolution(label,excitation_list=['rf1','rf2'])
        
        
        res_freq,res_signal,res_lockin,phase= lockin.lockin_detect(transient.data.conductivity[label][0:-2].size,
                                                             Simulation_parameters.time_intervale[label],
                                                             transient.data.conductivity[label][0:-2],
                                                             2)
        label='with resonance steady state fft'
        transient.data.current_operation_id+=1
        transient.data.operation_order[transient.data.current_operation_id]=label
        transient.data.time[label]=res_freq
        transient.data.conductivity[label]=res_signal
        transient.data.detected_lockin_signal.append([field_list[transition],res_lockin])
        
    #########################
    # Plot results
    print transient.data.detected_lockin_signal
    #transient.Plot_transient()# conductivity,light and resonance transient, fit/field
 
    clk.TOC("Total")
    clk.Display(level="L") #L/M/H
    
    
def Main_Lineshape(transition=0,scan_amplitude=0.002,scan_resolution=0.0001,name="scan"):
    #########################
    #Instanciation
    clk=Computation_time(True)
    clk.TIC("Total","H")
    data=Data("Report")
    sdr_param=SDR_parameters()
    Bi = Shallow_Donor(clk,"Bi",fortran_compiled=True)
    Si = Silicon(clk,fortran_compiled=True)
    lockin=Lockin_Detection(clk)
    
    field_list=[0.0660558,
                 0.0796147,
                 0.0986649,
                 0.12595,
                 0.164827,
                 0.217885,
                 0.285122,
                 0.363939,
                 0.450981,
                 0.543407]
    X=[]
    Y=[]
    Y_phase=[]
    
    rf1=ResonanceFrequency(clk,frequency=9076,power=1e10,modulated=False)
    rf2=ResonanceFrequency(clk,frequency=149.885,power=1e10,modulated=False)
    excitation={'rf1':rf1}
    
    b=Bfield(clk,field_list[transition],modulated=True)
    sdr_model=Parametrized_SDR_model(bulk=Si, donor=Bi,bfield=b, clk=clk,excitation_dic=excitation)
    transient=Transient(sdr_model=sdr_model,data=data,clk=clk)
    
    #########################
    #Initialization
    data.current_population=np.zeros([sdr_model.full_dim])
    data.current_population[sdr_model.c3i:sdr_model.c3e]=1/sdr_model.donor.dim_n*np.ones([sdr_model.donor.dim_n])

    #########################
    #Time evolution
    label='without resonance'
    transient.Temporal_evolution(label)
        
    for b0 in np.arange(field_list[transition]-scan_amplitude,field_list[transition]+scan_amplitude,scan_resolution):
        print "\n\nField= ",b0
        
        transient.sdr_model.bfield.B0=b0



        #########################
        #Time evolution
        #label='without resonance'
        #transient.Temporal_evolution(label)
        label='with resonance'
        transient.Temporal_evolution(label,excitation_list=['rf1'])
        label='with resonance steady state'
        transient.Temporal_evolution(label,excitation_list=['rf1'])
        
        
        res_freq,res_signal,res_lockin,res_lockin_phase= lockin.lockin_detect(transient.data.conductivity[label][0:-2].size,
                                                             Simulation_parameters.time_intervale[label],
                                                             transient.data.conductivity[label][0:-2],
                                                             2)
        label='with resonance steady state fft'
        transient.data.current_operation_id+=1
        transient.data.operation_order[transient.data.current_operation_id]=label
        transient.data.time[label]=res_freq
        transient.data.conductivity[label]=res_signal
        transient.data.detected_lockin_signal.append([b0,res_lockin])
        transient.data.detected_lockin_signal_phase.append([b0,res_lockin_phase])
        X.append(b0)
        Y.append(res_lockin)
        Y_phase.append(res_lockin_phase)
        
        
    #########################
    # Plot results
    print transient.data.detected_lockin_signal
    transient.data.Save_lockin(name,phase=True)
    
    fig = plt.figure()
    ax=fig.add_subplot(2,1,1)
    ax.plot(X,Y,'ob-')
    ax.set_title("2nd harmonic")
    ax.set_ylabel('SDR(a.u.)', color='b')
    for tl in ax.get_yticklabels():
        tl.set_color('b')
    
    ax2 = ax.twinx()
    ax2.plot(X, Y_phase, 'or-')
    ax2.set_ylabel('Phi', color='r')
    for tl in ax2.get_yticklabels():
        tl.set_color('r')
    
    plt.show()
        
    #transient.Plot_transient()# conductivity,light and resonance transient, fit/field
 
    clk.TOC("Total")
    clk.Display(level="L") #L/M/H
    
if __name__ == "__main__":
    #Simple transient
    #Main(transition=0,file_name="20121106_Bmod_cond7_2")
    Main_modfreq(transition=0,file_name="20121106_Bmod")
    #Main_Spectrum()
    
    '''
    Main_Lineshape(transition=0,
                   scan_amplitude=0.002,
                   scan_resolution=0.00005,
                   name="lineshape_20-1")
    '''

    